﻿using System;

namespace SoftwareConsulting.BI11.NumericalMethods.Interpolation
{    
    /// <summary>
    /// Класс интерполяции сплайнами
    /// </summary>
    /// <remarks>
    /// Может использоваться для интерполяции между рядами двумерных точек.
    /// Интерполирующая функция - кубический сплайн с существующими первыми производными в граничных точках.
    /// Если первые производные в граничных точках не существуют, то используется т.н. натуральный сплайн
    /// с нулевыми вторыми производными в граничных точках.
    /// Входные параметры:
    /// X   -   абсциссы точек;                    
    /// Y   -   значения функций в этих точках;                    
    /// N   -   число точек;
    /// T   -   точка, в которой осуществляется интерполяция
    /// Результат:
    /// значение интерполяционного полинома в точке T
    /// </remarks>
    public class SplineInterpolation 
    {
        private double _result;

        /// <summary>
        /// Конструктор по умолчанию
        /// </summary>                
        public SplineInterpolation(double[] x, double[] y, int N, double t)
        {
            _result = 0;
            double[] c = GetLinearSplineCoefficients(x, y, N);
            int n = 0;
            int l = 0;
            int r = 0;
            int m = 0;
            n = (int)Math.Round(c[2]);
            //Бинарный поиск в отрезке [ x[0], ..., x[n-2] ] (т.е. x[n-1] не включена)            
            l = 3;
            r = 3 + n - 2 + 1;
            while (l != r - 1)
            {
                m = (l + r) / 2;
                if (c[m] >= t)
                {
                    r = m;
                }
                else
                {
                    l = m;
                }
            }
            //собственно интерполяция
            t = t - c[l];
            m = 3 + n + 4 * (l - 3);
            _result = c[m] + t * (c[m + 1] + t * (c[m + 2] + t * c[m + 3]));
        }

        /*************************************************************************
   Построение таблицы коэффициентов кусочно-линейного сплайна

   Входные параметры:
       X           -   абсциссы, массив с нумерацией элементов [0..N-1]
       Y           -   значения функции,
                       массив с нумерацией элементов [0..N-1]
       N           -   число точек, N>=2

   Выходные параметры:
       C           -   таблица коэффициентов сплайна для использования в
                       подпрограмме SplineInterpolation
   *************************************************************************/
        /// <summary>
        /// Возвращает таблицу коэффициентов кусочно-линейного сплайна
        /// </summary>
        /// <param name="x">абсциссы точек</param>
        /// <param name="y">значения функций в этих точках</param>
        /// <param name="N">число точек</param>
        /// <returns>коэффициенты</returns>
        public double[] GetLinearSplineCoefficients(double[] x,double[] y, int N)
        {            
            int i = 0;
            int tblsize = 0;
            x = (double[])x.Clone();
            y = (double[])y.Clone();            
            SortPoints(ref x, ref y, N);
            tblsize = 3 + N + (N - 1) * 4;
            double[] c = new double[tblsize - 1 + 1];
            c[0] = tblsize;
            c[1] = 3;
            c[2] = N;
            for (i = 0; i <= N - 1; i++)
            {
                c[3 + i] = x[i];
            }
            for (i = 0; i <= N - 2; i++)
            {
                c[3 + N + 4 * i + 0] = y[i];
                c[3 + N + 4 * i + 1] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
                c[3 + N + 4 * i + 2] = 0;
                c[3 + N + 4 * i + 3] = 0;
            }
            return c;
        }

        /// <summary>
        /// Сортирует точки
        /// </summary>
        /// <param name="x">абсциссы точек</param>
        /// <param name="y">значения функций в этих точках</param>
        /// <param name="N">число точек</param>
        private void SortPoints(ref double[] x,ref double[] y,int n)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int t = 0;
            double tmp = 0;
            bool isascending = new bool();
            bool isdescending = new bool();
            isascending = true;
            isdescending = true;
            for (i = 1; i <= n - 1; i++)
            {
                isascending = isascending & x[i] > x[i - 1];
                isdescending = isdescending & x[i] < x[i - 1];
            }
            if (isascending)
                return;            
            if (isdescending)
            {
                for (i = 0; i <= n - 1; i++)
                {
                    j = n - 1 - i;
                    if (j <= i)
                        break;                    
                    tmp = x[i];
                    x[i] = x[j];
                    x[j] = tmp;
                    tmp = y[i];
                    y[i] = y[j];
                    y[j] = tmp;
                }
                return;
            }
            if (n == 1)
                return;            
            i = 2;
            do
            {
                t = i;
                while (t != 1)
                {
                    k = t / 2;
                    if (x[k - 1] >= x[t - 1])
                    {
                        t = 1;
                        continue;
                    }                    
                    tmp = x[k - 1];
                    x[k - 1] = x[t - 1];
                    x[t - 1] = tmp;
                    tmp = y[k - 1];
                    y[k - 1] = y[t - 1];
                    y[t - 1] = tmp;
                    t = k;                    
                }
                i = i + 1;
            }
            while (i <= n);
            i = n - 1;
            do
            {
                tmp = x[i];
                x[i] = x[0];
                x[0] = tmp;
                tmp = y[i];
                y[i] = y[0];
                y[0] = tmp;
                t = 1;
                while (t != 0)
                {
                    k = 2 * t;
                    if (k > i)
                    {
                        t = 0;
                        continue;
                    }                    
                    if (k < i)
                    {
                        if (x[k] > x[k - 1])
                        {
                            k = k + 1;
                        }
                    }
                    if (x[t - 1] >= x[k - 1])
                    {
                        t = 0;
                        continue;
                    }                    
                    tmp = x[k - 1];
                    x[k - 1] = x[t - 1];
                    x[t - 1] = tmp;
                    tmp = y[k - 1];
                    y[k - 1] = y[t - 1];
                    y[t - 1] = tmp;
                    t = k;                                        
                }
                i = i - 1;
            }
            while (i >= 1);
        }

        /// <summary>
        /// Возвращает решение
        /// </summary>        
        public double GetSolution()
        {
            return _result;
        }
    }
}
